auto_metrics Nightly Notebook¶
By Josh Dillon
Last Updated January 27, 2021
auto_metrics is a module in hera_qm that computes a series of statistics on night-long autocorrelation waterfalls to find antenna outliers in shape, power, or temporal structure. In general, these are assessed by collapsing each waterfall to a single (normalized) spectrum, then comparing that spectrum to the mean/median of the unflagged antennas' spectra to compute some difference metric. Those are then converted into modified Z-scores by comparing the overall distribution of good antennas and the worst antenna is flagged if it exceeds some threshold. This whole processes is repeated iteratively until no new bad antennas are identified. This proceeds in two rounds, first with more robust median-based statistics to identify the worst outliers, and then (after an RFI flagging step), a round with mean-based statistics. This notebook examines those mean-based spectra and statistics.
Statistics computed (after removing the worst offenders with a median-based metrics and then RFI flagging):¶
- Shape: Compute the mean spectrum over time for each autocorrelation, but then divide by the mean of the whole waterfall to get something near 1. Now compute the mean absolute value of the difference between that and the mean spectrum over all non-excluded antennas of the same polarization. Convert this to a modified Z-score by comparing to all non-excluded antennas.
- Power: Compute the mean spectrum over time for each autocorrelation without normalizing. Now compute the mean absolute value of the difference between the log of that (because we care about power outliers in dB rather than linear units) and the log of the mean spectrum over all non-excluded antennas of the same polarization. Convert this to a modified Z-score by comparing to all non-excluded antennas.
- Temporal variability: Divide each autocorrelation by the mean waterfall for all non-excluded antennas of the same polarization. Now reduce to a single spectrum by computing the standard deviation along the time axis for each. Now compute the mean value (not the absolute value, since low variability shouldn't get a high Z-score) of the difference between that and the mean spectrum over all non-excluded antennas of the same polarization. Convert this to a modified Z-score by comparing to all non-excluded antennas.
- Temporal discontinuities: Divide each autocorrelation by the mean waterfall for all non-excluded antennas of the same polarization. Now compute the element-by-element difference along the time axis, take the absolute value, and take the mean along the time axis to get a single spectrum. Now compute the mean value (not the absolute value, since low variability shouldn't get a high z-score) of the difference between that and the mean spectrum over all non-excluded antennas of the same polarization. Convert this to a modified Z-score by comparing to all non-excluded antennas.
import numpy as np
import pandas as pd
pd.set_option('display.max_rows', 1000)
import matplotlib.pyplot as plt
import matplotlib
import glob
import os
import operator
from hera_cal.io import HERAData
from hera_cal.utils import split_pol
from hera_qm.metrics_io import load_metric_file
from hera_notebook_templates.utils import status_colors
from IPython.display import display, HTML
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
Parse Inputs and Load Data¶
# If you want to run this notebook locally, copy the output of the next cell into the first few lines of this cell.
# JD = '2459122'
# data_path = '/lustre/aoc/projects/hera/H4C/2459122'
# os.environ["JULIANDATE"] = JD
# os.environ["DATA_PATH"] = data_path
# Use environment variables to figure out path to data
JD = os.environ['JULIANDATE']
data_path = os.environ['DATA_PATH']
print(f'JD = "{JD}"')
print(f'data_path = "{data_path}"')
print(f'Looking for data in {data_path} on JD {JD}...')
auto_metrics_file = sorted(glob.glob(os.path.join(data_path, f'zen.{JD}*.auto_metrics.h5')))
if len(auto_metrics_file) > 0:
auto_metrics_file = auto_metrics_file[0]
print(f'Found auto_metrics results file at {auto_metrics_file}.')
else:
raise OSError(f'{auto_metrics_file} not found.')
raw_auto_files = sorted(glob.glob(os.path.join(data_path, f'zen.{JD}.?????.sum.autos.uvh5')))
if len(raw_auto_files) > 0:
print(f'Found {len(raw_auto_files)} extracted autocorrelation files.')
else:
raise OSError(f'No files of the form zen.{JD}.?????.sum.autos.uvh5 found in {data_path}.')
# load auto_metrics and define use someful quantities
am = load_metric_file(auto_metrics_file)
mean_round_modz_cut = am['parameters']['mean_round_modz_cut']
ex_ants = am['ex_ants']['r2_ex_ants']
ants = sorted(set(bl[0] for bl in am['modzs']['r2_shape_modzs']))
# load raw autocorrelation waterfalls and define some useful quantities
hd = HERAData(raw_auto_files)
autos, _, _ = hd.read(axis='blt')
wf_shape = next(iter(autos.values())).shape
freqs = autos.freqs / 1e6
times = autos.times
lsts = autos.lsts * 12 / np.pi
# try to load a priori antenna statusesm but fail gracefully if this doesn't work.
a_priori_statuses = {ant: 'Not Found' for ant in ants}
try:
from hera_mc import cm_active
h = cm_active.ActiveData(at_date=JD)
h.load_apriori()
for ant_name in h.apriori:
ant = int("".join(filter(str.isdigit, ant_name)))
if ant in a_priori_statuses:
a_priori_statuses[ant] = h.apriori[ant_name].status
except Exception as err:
print(f'Could not load a priori antenna statuses.\nEncountered {type(err)} with message: {err}')
Summary Plots and Tables¶
# print ex_ants for easy copy-pasting to YAML file
print('ex_ants: [' + ", ".join(str(ant) for ant in ex_ants) + ']')
def Array_Plot():
plt.figure(figsize=(8, 8), dpi=75)
plt.scatter(np.array(list(autos.antpos.values()))[:,0],
np.array(list(autos.antpos.values()))[:,1], c='w', s=0)
for ant, pos in autos.antpos.items():
bad = ant in ex_ants
plt.gca().add_artist(plt.Circle(tuple(pos[0:2]), radius=7,
fill=(~bad), color=['grey','r'][bad]))
plt.text(pos[0],pos[1],str(ant), va='center', ha='center', color='w')
plt.xlabel("Antenna East-West Position (meters)")
plt.ylabel("Antenna North-South Position (meters)")
plt.title(f'Antenna Positions and Auto_Metrics Flags on {JD}\n(Maximum Modified Z-Score > {mean_round_modz_cut} in Red)');
plt.axis('equal')
plt.tight_layout()
Figure 1: Antenna Positions with auto_metrics flags.¶
This plot shows the antenna positions of all antennas in the data. The antennas with at least one Modified Z-score for one metric on one polarization exceeding the cut are entirely flagged.
Array_Plot()
# Parse modzs for Table 1 and other figures
modzs_to_check = {'Shape': 'r2_shape_modzs', 'Power': 'r2_power_modzs',
'Temporal Variability': 'r2_temp_var_modzs', 'Temporal Discontinuties': 'r2_temp_diff_modzs'}
worst_metrics = []
worst_zs = []
all_modzs = {}
binary_flags = {rationale: [] for rationale in modzs_to_check}
for ant in ants:
# parse modzs and figure out flag counts
modzs = {f'{pol} {rationale}': am['modzs'][dict_name][(ant, ant, pol)]
for rationale, dict_name in modzs_to_check.items() for pol in autos.pols()}
for pol in autos.pols():
for rationale, dict_name in modzs_to_check.items():
binary_flags[rationale].append(am['modzs'][dict_name][(ant, ant, pol)] > mean_round_modz_cut)
# figure out which metric is the largest outlier
worst_metric, worst_z = max(modzs.items(), key=operator.itemgetter(1))
worst_metrics.append(worst_metric)
worst_zs.append(worst_z)
# parse out all metrics for dataframe
for k in modzs:
col_label = k + ' Modified Z-Score'
if col_label in all_modzs:
all_modzs[col_label].append(modzs[k])
else:
all_modzs[col_label] = [modzs[k]]
# build dataframe
to_show = {'Ant': ants, 'A Priori Status': [a_priori_statuses[ant] for ant in ants],
'Worst Metric': worst_metrics, 'Worst Modified Z-Score': worst_zs}
to_show.update(all_modzs)
df = pd.DataFrame(to_show).sort_values('Worst Modified Z-Score', ascending=False)
# style dataframe
z_score_cols = [col for col in df.columns if col not in ['Ant', 'A Priori Status', 'Worst Metric']]
table1 = df.style.hide_index()\
.applymap(lambda val: 'font-weight: bold' if val in ex_ants else '', subset=['Ant']) \
.applymap(lambda val: 'color: red' if val in ex_ants else '', subset=['Ant']) \
.applymap(lambda val: f'background-color: {status_colors[val]}' if val in status_colors else '', subset=['A Priori Status']) \
.background_gradient(cmap='viridis', vmax=mean_round_modz_cut * 3, vmin=0, axis=None, subset=z_score_cols) \
.applymap(lambda val: 'font-weight: bold' if val > am['parameters']['mean_round_modz_cut'] else '', subset=z_score_cols) \
.applymap(lambda val: 'color: red' if val > am['parameters']['mean_round_modz_cut'] else '', subset=z_score_cols) \
.set_table_styles([dict(selector="th",props=[('max-width', '70pt')])])
Table 1: Modified Z-Score Summary¶
This table displays the metrics for each antenna, highlighting which one is the worst. It is sorted by each antenna's worst metric. When one metric exceeds the threshold, auto_metrics recommends cutting that antenna. Flagged antennas and metrics exceeding the cut are shown in bold and red. Also shown is the antenna's a priori status.
HTML(table1.render())
def Flag_Bar_Chart():
plt.figure(figsize=(8, 4), dpi=100)
# count
rationales = list(binary_flags.keys())
flags_list = np.array([binary_flags[rationale] for rationale in rationales])
antpol_flags = [np.sum(f) for f in flags_list]
ant_flags = [np.sum(np.array(f)[0::2] | np.array(f)[1::2]) for f in flags_list]
# make bar chart
plt.bar(np.arange(len(rationales)), antpol_flags, width=.7, color='yellow', ec='k', tick_label=rationales)
for x, (nflags, nants) in enumerate(zip(antpol_flags, ant_flags)):
plt.text(x, nflags/2, f'{nflags} Feeds on \n{nants} Antennas\nFlagged', va='center', ha='center')
# set labels
plt.ylabel('Antenna-Polarizations Flagged')
plt.xlabel('Reason for Flagging')
plt.tight_layout()
Figure 2: Flagging Rationale Summary¶
This bar chart summarizes the number of antenna-polarizations that are statistical outliers in each metric (though often they overlap). Some of these issues occur on both polarizations, so there are fewer unique antennas flagged for each rationale than there are ant-pols flagged, as noted by the labels.
Flag_Bar_Chart()
def Rationale_Corr_Plot():
plt.figure(figsize=(6,6), dpi=100)
# compute correlation matrix
rationales = list(binary_flags.keys())
flags_list = np.array([binary_flags[rationale] for rationale in rationales])
corrs = np.corrcoef(flags_list)
# plot and label correlation matrix
plt.imshow(corrs, cmap='viridis', interpolation='nearest', origin='upper')
for i in range(corrs.shape[0]):
for j in range(corrs.shape[1]):
plt.text(i, j, np.round(corrs[i, j], 3), va='center', ha='center',
bbox={'facecolor': 'w', 'ec': 'w', 'alpha': .75})
# colorbar, labels, and style
plt.yticks(range(len(rationales)), rationales)
plt.xticks(range(len(rationales)), rationales, rotation=-45, ha='right')
plt.gca().xaxis.tick_top()
plt.clim([0, 1])
plt.colorbar(fraction=0.046, pad=0.04)
plt.tight_layout()
Figure 3: Flagging Rationale Correlations¶
This plot shows the probability that if a given ant-pol is flagged for some reason, it's also flagged for another reason.
Rationale_Corr_Plot()
def plot_all_spectra(spectra, mod_zs, modz_cut, overall, freqs, reason, ex_ants=[],
xlabel='Frequency (MHz)', ylabel='', yscale='linear', ylim_factor=None):
'''Helper function for plotting all spectra and showing which ones were flagged and why.'''
fig, axes = plt.subplots(1,2, figsize=(14,5), dpi=100)
pols = sorted(set([bl[2] for bl in spectra]))
for ax, pol in zip(axes, pols):
# sort antennas into good, bad, and bad but not for this reason
bad_here = [bl for bl in spectra if (bl[2] == pol) and mod_zs[bl] > modz_cut]
other_bad = [bl for bl in spectra if (bl[2] == pol) and (bl[0] in ex_ants) and (bl not in bad_here)]
good = [bl for bl in spectra if (bl[2] == pol) and (bl[0] not in ex_ants)]
# plot all spectra
l1, l2, l3 = None, None, None
for bl in other_bad:
l2, = ax.plot(freqs, spectra[bl], 'darkviolet', lw=.5)
for bl in bad_here:
l1, = ax.plot(freqs, spectra[bl], 'r', lw=.5)
for bl in good:
l3, = ax.plot(freqs, spectra[bl], 'grey', alpha=.5)
l4, = ax.plot(freqs, overall[bl[2]], 'k--')
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_yscale(yscale)
if ylim_factor is not None:
ax.set_ylim([np.nanmin([spectra[bl] for bl in good]) / ylim_factor,
np.nanmax([spectra[bl] for bl in good]) * ylim_factor])
ax.set_title(f'Outliers in Autocorrelation {reason}: {pol}')
ax.legend([l1, l2, l3, l4], [f'Flagged for {reason}', 'Flagged for Another Reason', 'Unflagged', 'Average Unflagged'], loc='lower right')
plt.tight_layout()
overall_shape = {pol: np.nanmean([spec for bl, spec in am['spectra']['mean_spectra_normed'].items() if (bl[2] == pol)
and (bl[0] not in ex_ants)], axis=0) for pol in autos.pols()}
overall_power = {pol: np.nanmean([spec for bl, spec in am['spectra']['mean_spectra'].items() if (bl[2] == pol)
and (bl[0] not in ex_ants)], axis=0) for pol in autos.pols()}
overall_temp_var = {pol: np.nanmean([spec for bl, spec in am['spectra']['std_spectra_normed'].items() if (bl[2] == pol)
and (bl[0] not in ex_ants)], axis=0) for pol in autos.pols()}
overall_temp_diff = {pol: np.nanmean([spec for bl, spec in am['spectra']['mean_abs_diff_spectra_normed'].items() if (bl[2] == pol)
and (bl[0] not in ex_ants)], axis=0) for pol in autos.pols()}
Figure 4: Outliers in Autocorrelation Shape¶
This plot summarizes the spectra computed to compare to one another to find outliers in autocorrelation shape (see above for how that was computed). The mean compared to is shown as a black dashed line. Antennas in red were flagged as outliers, antennas in gray and purple were not. However, antennas in purple were flagged for some other reason, either another metric or on the other polarization. Completely flagged channels (RFI and band edges) appear as white gaps.
plot_all_spectra(am['spectra']['mean_spectra_normed'], am['modzs']['r2_shape_modzs'], mean_round_modz_cut,
overall_shape, freqs, 'Shape', ex_ants=ex_ants, yscale='linear', ylim_factor=1.2)
Figure 5: Outliers in Autocorrelation Power¶
This plot summarizes the spectra computed to compare to one another to find outliers in autocorrelation amplitude (see above for how that was computed). The mean compared to is shown as a black dashed line. Antennas in red were flagged as outliers, antennas in gray and purple were not. However, antennas in purple were flagged for some other reason, either another metric or on the other polarization. Completely flagged channels (RFI and band edges) appear as white gaps.
plot_all_spectra(am['spectra']['mean_spectra'], am['modzs']['r2_power_modzs'], mean_round_modz_cut,
overall_power, freqs, 'Power', ex_ants=ex_ants, yscale='log')
Figure 6: Outliers in Autocorrelation Temporal Variability¶
This plot summarizes the spectra computed to compare to one another to find outliers in autocorrelation temporal variability (as measured by a standard deviation over time; see above for how that was computed). The mean compared to is shown as a black dashed line. Antennas in red were flagged as outliers, antennas in gray and purple were not. However, antennas in purple were flagged for some other reason, either another metric or on the other polarization. Completely flagged channels (RFI and band edges) appear as white gaps.
plot_all_spectra(am['spectra']['std_spectra_normed'], am['modzs']['r2_temp_var_modzs'], mean_round_modz_cut,
overall_temp_var, freqs, 'Temporal Variability', ex_ants=ex_ants, yscale='log')
Figure 7: Outliers in Autocorrelation Temporal Discontinuities¶
This plot summarizes the spectra computed to compare to one another to find outliers in autocorrelation temporal discontinuities (as measured by the average absolute integration-to-integration difference over time; see above for how that was computed). The mean compared to is shown as a black dashed line. Antennas in red were flagged as outliers, antennas in gray and purple were not. However, antennas in purple were flagged for some other reason, either another metric or on the other polarization. Completely flagged channels (RFI and band edges) appear as white gaps.
plot_all_spectra(am['spectra']['mean_abs_diff_spectra_normed'], am['modzs']['r2_temp_diff_modzs'], mean_round_modz_cut,
overall_temp_diff, freqs, 'Temporal Discontinutities', ex_ants=ex_ants, yscale='log')
# compute average good autocorrelations for each polarization
avg_good_autos = {pol: np.zeros(wf_shape, dtype=float) for pol in autos.pols()}
for pol in autos.pols():
for i in range(wf_shape[0]):
avg_good_autos[pol][i] = np.mean([np.abs(autos[bl][i, :]) for bl in autos
if (bl[0] not in ex_ants) and (bl[2] == pol)], axis=0)
avg_good_autos[pol][am['flags']] = np.nan
def Avg_Auto_Plot():
fig, axes = plt.subplots(1, 2, figsize=(14,5), dpi=100)
for ax, pol in zip(axes, autos.pols()):
im = ax.imshow(avg_good_autos[pol], aspect='auto', interpolation='nearest',
extent=[freqs[0], freqs[-1], times[-1], times[0]])
ax.set_yticklabels(np.around(lsts[[min(max(np.searchsorted(times, t), 0), len(times) - 1) for t in ax.get_yticks()]], 2))
plt.colorbar(im, ax=ax)
ax.set_title(f'Average Good Raw {pol} Autocorrelation After Flagging')
ax.set_ylabel('LST (hours)')
ax.set_xlabel('Frequency (MHz)')
plt.tight_layout()
Figure 8: Average Good Autocorrelations and Flags¶
Here we show the waterfalls of the array-averaged autocorrelations over the night, after removing all flagged antennas. We also show the RFI mask generated between the median and mean rounds of antenna outlier detection. This is meant to show that there is little or no RFI remaining to affect the statistics.
Avg_Auto_Plot()
Per-Antenna Plots¶
def plot_spectra(axes, ant, spec, modzs, modz_cut, overall, reason, yscale='linear'):
'''Helper function for plotting both antennas of a given polarization.'''
for pol, ax in zip(autos.pols(), axes):
bl = (ant, ant, pol)
# plot good antennas
for bl2 in modzs:
if (bl2[0] not in ex_ants) and (bl2[2] == pol):
ax.plot(freqs, spec[bl2], 'grey', lw=.5, alpha=.5)
ax.plot(freqs, overall[pol], 'k--')
# plot this anetnna
color = 'r'
alpha = .75
if modzs[bl] >= modz_cut:
alpha = 1
elif bl[0] in ex_ants:
color = 'darkviolet'
else:
color = 'darkgreen'
ax.plot(freqs, spec[bl], color, alpha=alpha, label=f'{ant} {pol} (z = {np.round(modzs[bl],1)})')
# decorate axis
ax.set_yscale(yscale)
ax.legend(loc=4)
ax.set_title(f'{pol} {reason}')
len(ants)
def Plot_All_Auto_Metrics():
for row, ant in enumerate(df['Ant']):
# print title of section
display(HTML(f'<h2>Antenna {ant}:</h2>'))
# print metrics
df_row = df.loc[df['Ant'] == ant]
df_row = df_row.style.hide_index()\
.set_table_styles([dict(selector="th",props=[('max-width', '70pt')])])\
.applymap(lambda val: f'background-color: {status_colors[val]}' if val in status_colors else '', subset=['A Priori Status']) \
.applymap(lambda val: 'font-weight: bold' if val in ex_ants else '', subset=['Ant']) \
.applymap(lambda val: 'color: red' if val in ex_ants else '', subset=['Ant']) \
.applymap(lambda val: 'font-weight: bold' if val > am['parameters']['mean_round_modz_cut'] else '', subset=z_score_cols) \
.applymap(lambda val: 'color: red' if val > am['parameters']['mean_round_modz_cut'] else '', subset=z_score_cols)
display(HTML(df_row.render()))
# plot spectra and waterfalls
fig, axes = plt.subplots(2, 6, figsize=(18, 5.5), dpi=(50 * (104 / len(ants))**.5)) # this should help manage filesize
# plot individual spectra compared to all good antennas
plot_spectra(axes[:, 0], ant, am['spectra']['mean_spectra_normed'], am['modzs']['r2_shape_modzs'],
mean_round_modz_cut, overall_shape, 'Shape', yscale='linear')
plot_spectra(axes[:, 1], ant, am['spectra']['mean_spectra'], am['modzs']['r2_power_modzs'],
mean_round_modz_cut, overall_power, 'Power', yscale='log')
plot_spectra(axes[:, 2], ant, am['spectra']['std_spectra_normed'], am['modzs']['r2_temp_var_modzs'],
mean_round_modz_cut, overall_temp_var, 'Temporal Variability', yscale='linear')
plot_spectra(axes[:, 3], ant, am['spectra']['mean_abs_diff_spectra_normed'], am['modzs']['r2_temp_diff_modzs'],
mean_round_modz_cut, overall_temp_diff, 'Temporal Discontinutities', yscale='linear')
plt.tight_layout()
# plot linear-scale waterfalls
for pol, ax in zip(autos.pols(), axes[:, 4]):
bl = (ant, ant, pol)
im = ax.imshow(np.where(am['flags'], np.nan, autos[bl].real),
aspect='auto', interpolation='nearest',
extent=[freqs[0], freqs[-1], times[-1], times[0]])
ax.set_yticklabels(np.around(lsts[[min(max(np.searchsorted(times, t), 0), len(times) - 1) for t in ax.get_yticks()]], 2))
ax.set_title(f'{pol} Waterfall (Linear Scale)')
# plot log-scale mean-divided waterfalls
for pol, ax in zip(autos.pols(), axes[:, 5]):
bl = (ant, ant, pol)
to_plot = autos[bl].real / avg_good_autos[pol]
to_plot[am['flags']] = np.nan
to_plot /= np.nanmean(to_plot)
im = ax.imshow(np.log10(to_plot), aspect='auto', interpolation='nearest', vmin=-.07, vmax=.07,
extent=[freqs[0], freqs[-1], times[-1], times[0]])
ax.set_yticklabels(np.around(lsts[[min(max(np.searchsorted(times, t), 0), len(times) - 1) for t in ax.get_yticks()]], 2))
ax.set_title(f'{pol} Log(Normalized Waterfall)')
plt.show()
# print some whitespace
display(HTML('<hr style="height:3px">'))
Figure 9: Per-Antenna Statistics, Spectra, and Waterfalls¶
Here we show the metrics for each antenna and the spectra/waterfalls that hopefully explain what led to them. The table reproduces the information from Table 1 above. The first four panels in each row clearly highlight the antenna's spectrum as it compares to the mean good antenna (black) and the distribution of good antennas (gray). Spectra in red were flagged as outliers. Spectra in purple were flagged for some other reason, either another metric or on the other polarization. Good antennas are shown in green. Completely flagged channels (RFI and band edges) appear as white gaps. In the fifth column, the waterfall of that autocorrelation is shown on a linear scale after RFI/band edge flags (white). In the sixth column, we show the log (base 10) of the same waterfall, divided by the average good antennas' waterfall of that polarization and then normalized to an average of 1.
Plot_All_Auto_Metrics()
from hera_qm import version
print(version.construct_version_info())